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1.  Introduction 


This  report  documents  a  set  of  functions,  written  in  C++,  that  can  be  used  to  perform 
interpolations  (nearest- neighbor,  linear,  and  cubic)  and  to  find  coefficients  for  best-fit  equations. 
Functions  for  working  with  periodic  equations  are  included. 

Measurements  can  be  taken  and  stored  using  a  variety  of  schemes.  The  interpolations  described 
in  this  report  are  designed  to  work  with  data  sets  where  measured  values  for  all  independent  and 
dependent  variables  are  explicitly  stored.  Furthermore,  for  data  that  has  more  than  1  independent 
variable,  the  measurements  must  all  lie  on  some  type  of  rectangular  grid,  where  all  grid  locations 
are  populated. 

When  measurements  do  not  lie  on  some  type  of  rectangular  grid,  interpolations  become  more 
difficult.  For  those  types  of  data  sets,  KD-trees  can  be  used  to  perform  efficient  nearest-neighbor 
searches,  which  can  act  as  interpolations.  The  yKDTree  namespace^  can  be  used  to  work  with 
KD-trees  using  C++. 

The  yBilinear  namespace  can  be  used  for  the  special  case  of  data  sets  with  2  independent 
variables,  where  all  measurements  lie  on  a  rectangular  grid,  all  grid  locations  are  populated,  and 
all  grid  lines  are  evenly  spaced. 

The  functions  that  are  described  in  this  report  have  been  grouped  into  the  yinterp  namespace, 
which  is  summarized  at  the  end  of  this  report.  The  yinterp  namespace  relies  exclusively  on 
standard  C++  operations  and  functions.  However,  example  code  that  is  included  in  this  report 
makes  use  of  the  yRandom  namespace  for  generating  pseudorandom  numbers,  the  yBmp 
namespace"^  for  creating  images,  and  the  yIo2^  namespace  for  reading  and  parsing  text  files. 


2.  Background 


My  motivation  for  this  project  began  with  a  task  that  involved  interpolating  a  look-up  table  with 
8  independent  variables,  some  of  which  were  periodic.  Prior  to  this  project,  I  had  always  written 
interpolators  on  an  ad  hoc  basis.  However,  the  thought  of  having  to  deal  with  an  array  that 
requires  8  indices  to  access  a  particular  element  (something  like  A[i][j][k][l][m][n][o][p]) 
convinced  me  that  I  needed  a  more  formal  set  of  tools. 

The  task  required  that  all  code  be  written  in  C++.  To  maximize  portability,  the  code  had  to  be 
written  without  the  use  of  any  platform- specific  or  C++1 1  tools.  Performance  was  a  priority.  In 
addition,  I  wanted  a  versatile  set  of  tools  that  I  could  use  for  future  projects.  The  tools  needed  to 
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be  simple  to  use  and  well-tested  so  that  I  could  feel  comfortable  sharing  them  with  other 
programmers. 


3.  Searching  Sorted  Arrays:  The  Binary Search()  Function 


The  BinarySearchO  function  uses  a  binary  search  algorithm  to  search  arrays  that  are  sorted  into 
ascending  order.  For  arrays  that  do  not  contain  duplicate  entries,  the  BinarySearch()  function  can 
be  used  to  find  a  pointer  c  such  that 


*c  <  k  <  *(c-i-l),  (1) 

where  k  is  key  on  which  the  search  is  based. 

For  arrays  that  contain  duplicate  entries,  Eq.  1  needs  to  be  slightly  modified: 

*c<k<  *(c-i-l).  (2) 

Thus,  if  k  is  equal  to  an  element  that  has  one  or  more  duplicates,  c  may  point  to  any  one  of  the 
duplicates. 

Performance  for  the  BinarySearch()  function  is  0(log(n)),  which  is  demonstrated  in  Section  3.7. 

3.1  BinarySearchO  Code 

template<class  T>T*BinarySearch(//<======FIND  THE  POINTER  c  |  *c  <=  k  <  *(c+l) 

T*a,T*b,//<- -ARRAY  START  &  END  POINTS  (ARRAY  MUST  BE  SORTED  IN  INC.  ORDER) 


T  k){//< - KEY 

if (k<*a)return  a-1;// . note  that  a-1  may  point  to  an  invalid  address 


for(T*c; k<*--b; k>*c?a=c : b=c+l)c=a+(b-a)/2; /*->*/return  b; 

}//  YAGENAUT@GMAIL.COM - LAST~UPDATED~21DUL2014 - 


3.2  BinarySearchO  Template  Class 

T  T  must  be  a  sortable  data  type. 

3.3  BinarySearchO  Parameters 

a  a  points  to  the  beginning  of  the  array  that  will  be  searched.  The  array  should  be 

sorted  into  ascending  order. 

b  b  points  to  one  element  past  the  end  of  the  array  that  will  be  searched.  Thus,  b  is 

used  to  define  the  end  of  the  search  region  but  is  not  included  in  the  search,  b 
should  be  greater  than  a. 

k  k  is  the  key  on  which  the  search  is  based. 
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3.4  BinarySearchO  Return  Value 


For  arrays  that  do  not  contain  duplicate  entries,  the  BinarySearch()  function  returns  a  pointer  to 
the  greatest  array  element  that  is  less  than  or  equal  to  k.  If  k  is  less  than  the  least  array  element, 
then  a-1  is  returned.  Care  needs  to  be  taken  when  accessing  the  return  pointer,  since  a-1  may  not 
point  to  a  valid  address. 

3.5  BinarySearchO  Simple  Example 

The  following  example  uses  the  BinarySearchO  function  to  search  a  small  sorted  array  that 
contains  duplicates. 

#include  <cstdio>// . printf() 

#include  "y_interp. h"// . yinterp 

int  main(){//<=========A  SIMPLE  EXAMPLE  FOR  THE  yinterp: :BinarySearch( )  FUNCTION 

int  X[ll] ={30, 40, 41, 42, 42, 42, 42, 43, 44, 45, 50}; 

printf("X={");/*&*/for(int  i=0; i<ll;++i)printf ( "%d%s",X[i] , i ! =10?", " : "}\n\n" ) ; 

printfC'  *c  I  KEY  |  *(c+l)  |  INDEX\n - \n"); 

for(int  k=28; k<52;++k){ 

int*c=ylnterp: : BinarySearch(X,X+ll,  k) ; 
if(c<X)printf("  --  |%4d  |%5d  |%4d\n",k,*(c+l),c-X); 

else  if (c==X+10)printf ( "%3d  |%4d  |  --  |%4d\n",*c,k,c-X); 

else  printf("%3d  |%4d  |%5d  | %4d\n", *c, k, *(c+l) , c-X) ; } 

^11  E  N  AUT  ^^GI^A  XL*  LAS  T^U  PDAT  ED^2XDLJL20 


OUTPUT: 


X={30 

,40,41 

,42,42,42,42,43,44,45,50} 

*c  1 

KEY  1 

*(c+l) 

1  INDEX 

1 

28  1 

30 

1  -1 

--  1 

29  1 

30 

1  -1 

30  1 

30  1 

40 

1  3 

30  1 

31  1 

40 

1  3 

30  1 

32  1 

40 

1  3 

30  1 

33  1 

40 

1  3 

30  1 

34  1 

40 

1  3 

30  1 

35  1 

40 

1  3 

30  1 

36  1 

40 

1  3 

30  1 

37  1 

40 

1  3 

30  1 

38  1 

40 

1  3 

30  1 

39  1 

40 

1  3 

40  1 

40  1 

41 

1  1 

41  1 

41  1 

42 

1  2 

42  1 

42  1 

42 

1  5 

43  1 

43  1 

44 

1  7 

44  1 

44  1 

45 

1  3 

45  1 

45  1 

50 

1  9 

45  1 

46  1 

50 

1  9 

45  1 

47  1 

50 

1  9 

45  1 

48  1 

50 

1  9 

45  1 

49  1 

50 

1  9 

50  1 

50  1 

-- 

1  10 

50  1 

51  1 

-- 

1  10 
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3.6  BinarySearchO  Text  Example 

The  following  example  uses  the  BinarySearch()  function  to  determine  where  a  word  should  be 
inserted  into  a  sorted  list. 


#include  <iostream>// . cout 

#include  <string>// . string 

#include  "y_interp. h"// . yinterp 

int  main(){//<===========A  TEXT  EXAMPLE  FOR  THE  yinterp: :BinarySearch()  FUNCTION 

std : : string  S[6]={ "defective", "defend" , "defendant", "defender", "defense" , 
"defensive" }, s="defenest rate" ; 
for(int  i=0; i<6;++i)std : : cout<<S[i]<< (i==5? "\n\n" : ",  "); 
std: : string*c=ylnterp: :BinarySearch(S,S+6,s); 

std: :cout<<"The  word  "<<s<<"  should  be  between  "<<*c<<"  and  "<<*(c+l)<<" . \n"; 
^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  PDAT  EO'^21DUL20 


OUTPUT: _ 

defective,  defend,  defendant,  defender,  defense,  defensive 

The  word  defenestrate  should  be  between  defender  and  defense. 


3.7  BinarySearchO  Performance 

The  following  example  begins  by  using  the  yRandom  namespace  to  populate  an  array  with  2^"^ 
pseudorandom  numbers  using  the  Mersenne  twister  19937  algorithm.  The  array  is  then  sorted 
using  the  sort()  function.  Average  search  times  are  calculated  for  searches  performed  using  a 
linear  search  algorithm  and  the  BinarySearchO  function.  Figure  1  presents  a  graph  of  the  results. 

The  average  search  time  for  the  LinearSearchO  function  is  0(n) .  The  average  search  time  for  the 
BinarySearchO  function  is  0(lo^n)). 
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#include  <algorithm>// . sort() 

#include  <cstdio>// . printf() 

#include  <ctime>// . clock( ) ,CLOCKS_PER_SEC 

#include  "y_interp. h"// . yinterp 

#include  "y_random. h"// . yRandom 

template<class  T>T*LinearSearch(//<========FIND  THE  POINTER  c  |  *c  <=  k  <  *(c+l) 

T*a,T*b,//< - ARRAY  START  &  END  POINTS  (ARRAY  MUST  BE  SORTED  IN  INC.  ORDER) 

T  k){//< - KEY 

if (k<*a)return  a-1;// . note  that  a-1  may  point  to  an  invalid  address 

while(k<*--b);/*->*/return  b; 


}// 


YAGENAUT@GMAIL.COM 


LAST~UPDATED~21DUL2014 


int  main(){//<=========PERFORMANCE  TEST  FOR  THE  yinterp: :BinarySearch()  FUNCTION 

const  int  N=16384,M=10000000; // . . .max  #  of  elements  in  array,  #  of  repetitions 
unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 

double*X=new  double[N];/*<-*/for(int  i=0;i<N;++i)X[i]=yRandom: :RandU(I,0,N); 
std : : sort(X,X+N) ; // . X  is  a  sorted,  unchanging  list  of  numbers  to  search 


linear  search 


binary 


size 

of 

array 


search 
time  (s) 


index 

average 


search 
time  (s) 


search\n" ) ; 

- \n"); 

index\n" ) ; 
average\n" ) ; 
- \n"); 


printf (' 
printf (' 
printf (' 
printf (' 
printf (' 
for(int  m=2;m<16385;m*=2){ 

printf("%8d  | ",m), yRandom: :Initialize(I,l); 
double  y=0, z=0,t=clock( ) ; 

for (int  i=0; i<M;++i)y+=LinearSearch(X,X+m, yRandom: : RandU(I,0,m) ) -X; 
printf("%8.3f  |%12.6f  | " , (clock( ) -t)/CLOCKS_PER_SEC,y/M) ,t=clock( ) ; 
yRandom: :Initialize(I,l); 

for(int  i=0;i<M;++i)z+=yInterp: : BinarySearch(X,X+m, yRandom: :RandU(I,0,m)) -X; 
printf ( "%8 . 3f  | %12 . 6f \n" , (clock( ) -t )/CLOCKS_PER_SEC, z/M) ; } 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  PDAT  E  1 D  U  L  2  0 


OUTPUT: 


size 

of 

linear  search  | 

1 

binary 

search 

search 

1 

1  index  | 

search  | 

index 

array 

time  (s) 

1  average  | 

1 _ 1 

time  (s)  1 

average 

2 

0.140 

1  -0.731079  1 

0.125  1 

-0.731079 

4 

0.140 

1  0.555495  1 

0.156  1 

0.555495 

8 

0.203 

1  1.833760  1 

0.224  1 

1.833760 

16 

0.244 

1  5.133788  1 

0.234  1 

5.133788 

32 

0.312 

1  13.015384  1 

0.327  1 

13.015384 

64 

0.390 

1  31.287311  1 

0.359  1 

31.287311 

128 

0.531 

1  66.159459  j 

0.421  1 

66.159459 

256 

0.936 

1  128.260683  j 

0.530  1 

128.260683 

512 

1.545 

1  258.242078  j 

0.592  1 

258.242078 

1024 

2.902 

1  519.408822  j 

0.671  1 

519.408822 

2048 

5.678 

1  1035.587832  j 

0.733  1 

1035.587832 

4096 

11.217 

1  2054.054591  j 

0.795  1 

2054.054591 

8192 

22.142 

1  4091.340401  1 

0.873  1 

4091.340401 

16384 

43.945 

1  8193.415153  j 

0.968  1 

8193.415153 
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Fig.  1  Performance  for  linear  and  binary  searches  (right  graph  is  semi-log) 


4.  Searching  Sorted  Arrays:  The  PeriodicSearchQ  Function 


The  PeriodicSearchO  function  uses  the  BinarySearch()  function  to  search  arrays  that  are 
associated  with  periodic  functions,  such  as  /  ,  where 

fix  +  np)  =  fix) .  (3) 

Thus,  the  PeriodicSearchO  function  can  be  used  to  find  a  pointer  c  such  that 

*c  <  k -I- np  < *  *(c-i-l),  (4) 

where  k  and  p  are  user  supplied  values  and  the  integer  n  is  chosen  such  that 

*a  <  k  -I-  np  <  *a-i-p.  (5) 


4.1  PeriodicSearchO  Code 

template<class  T>T*PeriodicSearch(//<=======FOR  ARRAYS  WITH  PERIODIC  VARIABLES 

T*a,T*b,//<- -STARTING  &  ENDING  POINTS  (ARRAY  MUST  BE  SORTED  IN  INC.  ORDER) 

T&k,//< - KEY  (WILL  BE  SET  TO  k'=k+np) 

T  p){//< - PERIOD  (p>0) 

return  BinarySearch(a, b, k=fmod(k-*a, p)+*a+(k-*a<0?p: 0) ) ; 

}//  YAGENAUT@GMAIL.COM - LAST~UPDATED~21DUL2014 - 


4.2  PeriodicSearchO  Template  Class 

T  T  must  be  a  sortable  data  type. 

4.3  PeriodicSearchO  Parameters 

a  a  points  to  the  beginning  of  the  array  that  will  be  searched.  The  array  should  be 

sorted  into  ascending  order. 
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b  b  points  to  one  element  past  the  end  of  the  array  that  will  be  searched.  Thus,  b  is 

used  to  define  the  end  of  the  search  region  but  is  not  included  in  the  search,  b 
should  be  greater  than  a. 

k  k  is  the  key  on  which  the  search  is  based.  Note  that  before  the  function  returns,  k 

is  set  to  k',  where  k'  =  k  +  np. 

p  p  is  the  period  of  the  function  associated  with  the  array.  Values  for  p  are  typically 

chosen  such  that  p  >  *(b-l)  -  *a.  In  no  case  should  p  be  equal  to  zero. 

4.4  PeriodicSearchO  Return  Value 

The  PeriodicSearchO  function  returns  a  pointer  that  satisfies  the  requirements  for  c  in  Eq.  4. 

4.5  PeriodicSearchO  Example 

The  following  example  begins  by  creating  a  pair  of  arrays  that  are  used  to  store  independent  and 
dependent  variables  that  approximate  Eq.  6. 

y{x)-  sin(2;w  /  5)  + 1 .  (6) 

The  BinarySearchO  and  PeriodicSearchO  functions  are  then  used  to  retrieve  dependent  variables 
based  on  a  selection  of  independent  variables.  Eigure  2  presents  a  pair  of  graphs  that  display  the 
results  of  the  searches. 


#include  <cstdio>// . f reopen ( ) ,  printf ( ) ,  stdout 

#include  "y_interp. h"// . yinterp, <cmath>{sin( )} 

int  main( ){//<==============COMPARISON  BETWEEN  BinarySearch( )  &  PeriodicSearchO 

f reopen ( "sine. csv'O  "w", stdout) ; 

double  X[20] ; /*<-*/for(int  i=0;i<20;++i)X[i]=5*i/20.+5; 

double  Y[20] ; /*<-*/for(int  i=0; i<20;++i)Y[i]=sin(2*3 . 14159*X[i]/5)+l; 

printf("x,y  -  BinarySearch( ) ),y  -  PeriodicSearch()\n"); 

for(double  x=-5,k;x<15;x+=.l){ 

printf ( 'Ox,Y[yInterp: :BinarySearch(X+l,X+20,x) -X] ); 
printf ( "%f\n",Y[yInterp: : PeriodicSearch(X,X+20, k=x, 5 . ) -X] ) ; } 

}//—  rf/vi/N/YAGENAUT@GMAI  L . 


7 


BinarySearchO 

... 

— 

5 

3 

5  1 

1  PeriodicSearchO 

0  I 

^  A 

5 

■■ . ■■■ 

... 

■■ 

■n  ■■ 

" 

5 

)  5  10  15 

Fig.  2  A  comparison  of  the  Search()  and  PeriodicSearch()  functions 


5.  Performing  Nearest-Neighbor  Interpolations:  The  NNInterpO  Function 


Suppose  that  some  function  y(x)  is  represented  by  a  finite  set  of  points  (x, ,  y, ),  such  as  is 
shown  in  Fig.  3. 


Fig.  3  y(x)  represented  by  a  finite  number  of  points,  with  nearest-neighbor  interpolation  shown  in  red 

The  NNInterpO  function  uses  Eq.  7,  represented  by  the  solid  red  line  in  Fig.  3,  to  approximate 

yW: 
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y(x)  =  < 


(V) 


y,+i  forx>^ 


y.  otherwise 


5.1  NNInterpO  Code 

template<class  T>T  NNInterp(//<==================NEAREST-l\IEIGHBOR  INTERPOLATOR 


const  T*X,//< - BRACKETING  X  VALUES  (*X  AND  X[l]  MUST  BE  VALID) 

const  T*Y,//< - BRACKETING  Y  VALUES  (*Y  AND  Y[l]  MUST  BE  VALID) 

T  x){//< - VALUE  TO  INTERPOLATE  AT  (TYPICALLY,  *X  <=  x  <=  x[l]) 


return  x>(*X+X[l] )/2?Y[l] : *Y; 

}//  YAGENAUT@GMAIL.COM - LAST~UPDATED~21DUL2014 - 


5.2  NNInterpO  Template  Class 

T  T  is  typically  a  floating-point  data  type. 

5.3  NNInterpO  Parameters 

X  X  points  to  an  array  that  is  used  to  store  values  for  x,  from  Eq.  7.  Specifically, 

*X  =  X-  and  X[l]  =  .  Thus,  both  X  and  X-i-1  must  point  to  valid  addresses. 

Y  Y  points  to  an  array  that  is  used  to  store  values  for  y.  from  Eq.  7.  Specifically, 

*Y  =  y^  and  Y[l]  =  .  Thus,  both  Y  and  Y-i-1  must  point  to  valid  addresses. 

X  X  represents  the  independent  variable  x  from  Eq.  7. 

5.4  NNInterpO  Return  Value 

The  NNInterpO  function  returns  y  from  Eq.  7. 

5.5  NNInterpO  Simple  Example 

The  following  example  begins  by  creating  a  pair  of  arrays  that  are  used  to  store  independent  and 
dependent  variables  that  approximate  Eq.  6.  Next,  the  BinarySearchO  function  is  used  to  find  a 
pointer  that  conforms  to  Eq.  1,  which  is  then  converted  to  an  index.  Einally,  the  NNInterpO 
function  is  used  to  approximate  y(7.I80) . 

Note  that  the  value  for  the  BinarySearchO  parameter  a  is  set  to  X-i-1  to  avoid  the  possibility  of 
accessing  the  X-1  element.  Similarly,  the  value  for  the  BinarySearchO  parameter  b  is  set  to 
X-I-1 9  to  avoid  the  possibility  of  accessing  the  X-i-20  element. 
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#include  <cstdio>// . printf() 

#include  "y_interp. h"// . yinterp, <cmath>{sin( )} 

int  main(){//<=============A  SIMPLE  EXAMPLE  FOR  THE  yinterp: :NNInterp( )  FUNCTION 

double  X[20] ; /*<-*/for(int  i=0;i<20;++i)X[i]=5*i/20.+5; 

double  Y[20] ; /*<-*/for(int  i=0; i<20;++i)Y[i]=sin(2*3 . 14159*X[i]/5)+l; 

double  x=7.18; 

int  i=ylnterp: :BinarySearch(X+l,X+19,x)-X; 

double  y=ylnterp: :NNInterp(X+i,Y+i,x); 

printfC'At  x=%.3f,  y  is  approximately  %.3f . \n",x,y); 

"^11  E  N  AUT  ^)GI^AI  L  •  COI^*^*^^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^  LAS  T^U  PD  AT  E  D^2  T  D  U  L  2  0 


OUTPUT: _ 

At  x=7.180,  y  is  approximately  1.309. 


6.  Performing  Linear  Interpolations:  The  LinInterpO  Function 


Suppose  that  some  function  }:(x)  is  represented  by  a  finite  set  of  points  (x,  such  as  is 
shown  in  Fig.  4. 


Fig.  4  y{x)  represented  by  a  finite  number  of  points,  with  linear  interpolation  shown  in  red 

The  LinInterpO  function  uses  Eq.  8,  represented  by  the  solid  red  line  in  Fig.  4,  to  approximate 
y(x): 
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X  —  Xj 


(8) 


yix)  =  y,  +(y,^i 


X., 


-X, 


6.1  LinInterpO  Code 

template<class  T>T  Linlnterp(//<===========================LINEAR  INTERPOLATOR 


const  T*X,//< - BRACKETING  X  VALUES  (*X  AND  X[l]  MUST  BE  VALID) 

const  T*Y,//< - BRACKETING  Y  VALUES  (*Y  AND  Y[l]  MUST  BE  VALID) 

T  x){//< - VALUE  TO  INTERPOLATE  AT  (TYPICALLY,  *X  <=  x  <=  x[l]) 


return*Y+(Y[l]-*Y)*(x-*X)/(X[l]-*X); 

}//  YAGENAUT@GMAIL.COM - LAST~UPDATED~21DUL2014 - 


6.2  LinInterpO  Template  Class 

T  T  is  typically  a  floating-point  data  type. 

6.3  LinInterpO  Parameters 

X  X  points  to  an  array  that  is  used  to  store  values  for  x,  from  Eq.  8.  Specifically, 

*X  =  X,  and  X[l]  =  x,^j .  Thus,  both  X  and  X-i-1  must  point  to  valid  addresses. 

Y  Y  points  to  an  array  that  is  used  to  store  values  for  y.  from  Eq.  8.  Specifically, 

*Y  =  y.  and  Y[l]  =  .  Thus,  both  Y  and  Y-i-1  must  point  to  valid  addresses. 

X  X  represents  the  independent  variable  x  from  Eq.  8. 

6.4  LinInterpO  Return  Value 

The  LinInterpO  function  returns  y  from  Eq.  8. 

6.5  LinInterpO  Simple  Example 

The  following  example  begins  by  creating  a  pair  of  arrays  that  are  used  to  store  independent  and 
dependent  variables  that  approximate  Eq.  6.  Next,  the  BinarySearchO  function  is  used  to  find  a 
pointer  that  conforms  to  Eq.  1,  which  it  then  converted  to  an  index.  Einally,  the  LinInterpO 
function  is  used  to  approximate  y(7.I80) . 

Note  that  the  value  for  the  BinarySearchO  parameter  a  is  set  to  X-i-1  to  avoid  the  possibility  of 
accessing  the  X-1  element.  Similarly,  the  value  for  the  BinarySearchO  parameter  b  is  set  to 
X-I-1 9  to  avoid  the  possibility  of  accessing  the  X-i-20  element. 
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#include  <cstdio>// . printf() 

#include  "y_interp. h"// . yinterp, <cmath>{sin( )} 

int  main(){//<============A  SIMPLE  EXAMPLE  FOR  THE  yinterp: : Linlnterp( )  FUNCTION 

double  X[20] ; /*<-*/for(int  i=0;i<20;++i)X[i]=5*i/20.+5; 

double  Y[20] ; /*<-*/for(int  i=0; i<20;++i)Y[i]=sin(2*3 . 14159*X[i]/5)+l; 

double  x=7.18; 

int  i=ylnterp: :BinarySearch(X+l,X+19,x)-X; 

double  y=ylnterp: : LinInterp(X+i,Y+i,x); 

printfC'At  x=%.3f,  y  is  approximately  %.3f . \n",x,y); 

"^11  E  N  AUT  ^)GI^AI  L  •  COI^*^*^^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^  LAS  T^U  PD  AT  E  D^2  T  D  U  L  2  0 


OUTPUT: _ 

At  x=7.180,  y  is  approximately  1.387. 


7.  Performing  Cubic  Interpolations:  The  CubeInterpO  Function 


A  cubic  Hermite  spline  is  a  third-degree-polynomial  interpolating  function  that  is  uniquely 
determined  by  2  endpoint  positions  ( /ig  and  Py )  and  tangent  vectors  at  the  2  endpoint  positions 
(mg  and  mj ,  respectively): 

//g(t) ,  H^{t) ,  H^it) ,  and  are  known  as  Hermite  basis  functions  and  are  given  by  Eqs. 
10-13.  t  is  the  linearly  scaled  distance  from  p^  to  p^  such  that  at  p^,  t  =  0  and  at 

Pi,t  =  l. 


//g(0  =  2t^-3t"+l. 

(10) 

(11) 

H^{t)  =  -2t^  +  3t^ . 

(12) 

(13) 

Suppose  that  some  function  y(x)  is  represented  by  a  finite  set  of  points  (x,  such  as  is 
shown  in  Fig.  5. 
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Fig.  5  y(x)  represented  by  a  finite  number  of  points,  with  cubic  interpolation  shown  in  red 


From  Eq.  9  (and  noting  that  =  1  - H2(t) ), 


at 


at 


yf 


(14) 


Eq.  15  can  be  used  to  eonvert  from  the  sealed  independent  variable  t  to  the  general  independent 
variable  x . 


t  = 


x-x, 


(15) 


Eq.  16,  along  with  the  chain  rule,  ean  be  used  to  convert  the  derivatives. 

dt  _  1 

dx  x,^.j  -  x^ 


Thus, 


and 


(x,.^i  -x,.)m,.. 


dy 


dt 


1=1 


(x,.^i  -x,)m,.^i. 


(16) 


(17) 


(18) 
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where 


m..  = 


dy 

dx 


(19) 


Substituting  Eqs.  11-13,  17,  and  18  into  Eq.  14,  then  regrouping  terms,  results  in  the  form  of  the 
eubie  interpolating  function  that  is  used  by  the  CubeInterpO  function: 

y  =  y,.  +  (x  -  x,.)[{(m.  +  -  (2m,  +  +  m,]  +  (y,.^j  -  y,.)(3  -  2t)t^ .  (20) 

7.1  CubeInterpO  Code 

template<class  T>T  Cubelnterp(//<==========CUBIC  (HERMITE  SPLINE)  INTERPOLATOR 


const  T*X,//< - BRACKETING  X  VALUES  (*X  AND  X[l]  MUST  BE  VALID) 

const  T*Y,//< - BRACKETING  Y  VALUES  (*Y  AND  Y[l]  MUST  BE  VALID) 

T  x,//< - VALUE  TO  INTERPOLATE  AT  (TYPICALLY,  *X  <=  x  <=  x[l]) 

T  m0,T  ml){//< - SLOPES  AT  *X  AND  X[l] 


T  t=(x-*X)/(X[l]-*X); 

return*Y+(x-*X)*( ( (m0+ml)*t- (2*m0+ml) )*t+m0)+(Y[l] -*Y)*(3-2*t)*t*t; 
}//  YAGENAUT@GMAI  L  . 


7.2  CubeInterpO  Template  Class 

T  T  is  typically  a  floating-point  data  type. 

7.3  CubeInterpO  Parameters 

X  X  points  to  an  array  that  is  used  to  store  values  for  x,  from  Eq.  20.  Specifically, 

*X  =  X,  and  X[l]  =  x,^j .  Thus,  both  X  and  X-i-1  must  point  to  valid  addresses. 

Y  Y  points  to  an  array  that  is  used  to  store  values  for  y.  from  Eq.  20.  Specifically, 

*Y  =  y.  and  Y[l]  =  y,.^j .  Thus,  both  Y  and  Y-i-1  must  point  to  valid  addresses. 

mO  mO  represents  m, ,  the  slope  of  the  interpolating  function  at  x, . 

ml  ml  represents  ,  the  slope  of  the  interpolating  function  at  x,^j . 

X  X  represents  the  independent  variable  x  from  Eq.  20. 

7.4  CubeInterpO  Return  Value 

The  CubeInterpO  function  returns  y  from  Eq.  20. 

7.5  CubeInterpO  Simple  Example 

The  following  example  begins  by  creating  a  pair  of  arrays  that  are  used  to  store  independent  and 
dependent  variables  that  approximate  Eq.  6.  Next,  the  BinarySearch()  function  is  used  to  find  a 
pointer  that  conforms  to  Eq.  1,  which  is  then  converted  to  an  index.  Einally,  the  CubeInterpO 
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function  is  used  to  approximate  y(7.180) .  For  simplicity,  m,  and  have  been  set  to  zero. 
Section  8  presents  a  function  that  can  be  used  to  find  values  for  and  . 

Note  that  the  value  for  the  BinarySearch()  parameter  a  is  set  to  X+1  to  avoid  the  possibility  of 
accessing  the  X-1  element.  Similarly,  the  value  for  the  BinarySearch()  parameter  b  is  set  to 
X+1 9  to  avoid  the  possibility  of  accessing  the  X+20  element. 


#include  <cstdio>// . printf() 

#include  "y_interp. h"// . yinterp, <cmath>{sin( )} 

int  main(){//<===========A  SIMPLE  EXAMPLE  FOR  THE  yinterp; :CubeInterp()  FUNCTION 

double  X[20] ; /*<-*/for(int  i=0;i<20;++i)X[i]=5*i/20.+5; 

double  Y[20] ; /*<-*/for(int  i=0; i<20;++i)Y[i]=sin(2*3 . 14159*X[i]/5)+l; 

double  x=7.18; 

int  i=ylnterp: :BinarySearch(X+l,X+19,x)-X; 
double  y=ylnterp: :CubeInterp(X+i,Y+i,x,0. ,0. ); 
printfC'At  x=%.3f,  y  is  approximately  %.3f . \n",x,y); 

^11  YAG  E  N  AUT  ^^GI^AI  L  •  LAS  T  PDAT  EO'^21DUL20 


OUTPUT: 

At  x=7.180,  y  is  approximately  1.362. 


8.  Cardinal  Splines:  The  CardinalSlopeO  Function 


The  CardinalSlopeO  function  uses  Eq.  21  to  calculate  the  CubeInterp()-function  input  parameters 

mO  and  ml. 


m,  =(l-0^'^  ■  (21) 

-^,+1 

t  is  sometimes  referred  to  as  a  tension  parameter  and  is  typically  limited  to  the  interval  [0,1]. 

8.1  CardinalSlopeO  Code 

template<class  T>T  CardinalSlope(//<========CALCULATES  SLOPES  FOR  Cubelnterp() 


const  T*X,//< - BRACKETING  X  VALUES  (X[-l]  AND  X[l]  MUST  BE  VALID) 

const  T*Y,//< - BRACKETING  Y  VALUES  (Y[-l]  AND  Y[l]  MUST  BE  VALID) 

T  t){//< - TENSION  PARAMETER 


return(l-t)*(Y[l]-Y[-l])/(X[l]-X[-l]); 

}//  YAGENAUT@GMAIL.COM - LAST~UPDATED~213UL2014 - 


8.2  CardinalSlopeO  Template  Class 

T  T  is  typically  a  floating-point  data  type. 
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8.3  CardinalSlopeO  Parameters 

X  X  points  to  an  array  that  is  used  to  store  values  for  x,  from  Eq.  21.  Speeifioally, 

X[-l]  =  x-_^  and  X[l]  =  x-^j .  Thus,  both  X-1  and  X+1  must  point  to  valid 
addresses. 

Y  Y  points  to  an  array  that  is  used  to  store  values  for  y.  from  Eq.  21.  Speeifieally, 

Y[-l]  =  y,_i  and  Y[l]  =  y._^^ .  Thus,  both  Y-1  and  Y+1  must  point  to  valid 
addresses. 

t  t  represents  t ,  the  tension  parameter  from  Eq.  21 . 

Beeause  the  CardinalSlopeO  funetion  looks  at  points  to  either  side  of  the  braeketing  points,  extra 
care  needs  to  be  taken  when  selecting  X  and  Y. 

8.4  CardinalSlopeO  Return  Value 

The  CardinalSlopeO  function  returns  m.  from  Eq.  21. 

8.5  CardinalSlopeO  Simple  Example 

The  following  example  begins  by  creating  a  pair  of  arrays  that  are  used  to  store  independent  and 
dependent  variables  that  approximate  Eq.  6.  Next,  the  BinarySearchO  function  is  used  to  find  a 
pointer  that  conforms  to  Eq.  1,  which  is  then  converted  to  an  index.  Einally,  the  CubeInterpO 
function  is  used  to  approximate  y(7.180) . 

Note  that  the  value  for  the  BinarySearchO  parameter  a  has  been  set  to  X+1  to  avoid  the 
possibility  of  accessing  the  X-1  element.  Similarly,  the  value  for  the  BinarySearchO  parameter  b 
has  been  set  to  X+19  to  avoid  the  possibility  of  accessing  the  X+20  element. 

#include  <cstdio>// . printf() 

#include  "y_interp. h"// . ylnterp,<cmath>{sin( ) } 

int  main(){//<========A  SIMPLE  EXAMPLE  FOR  THE  yinterp: :CardinalSlope( )  FUNCTION 

double  X[20] ; /*<-*/for(int  i=0;i<20;++i)X[i]=5*i/20.+5; 

double  Y[20] ; /*<-*/for(int  i=0; i<20;++i)Y[i]=sin(2*3 . 14159*X[i]/5)+l; 

double  x=7.18; 

int  i=ylnterp: :BinarySearch(X+l,X+19,x)-X; 
double  m0=i==0| | i>17?0:ylnterp: :CardinalSlope(X+i,Y+i,0. ), 
ml=i==0 1  I i >17 ?0: yinterp: :CardinalSlope(X+i+l, Y+i+1,0. ) ; 
double  y=ylnterp: :CubeInterp(X+i,Y+i,x,m0,ml); 
printfC'At  x=%.3f,  y  is  approximately  %.3f .  \n'0x,y); 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  PDAT  ED'^21DUL20 


OUTPUT: _ 

At  x=7.180,  y  is  approximately  1.391. 
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9.  Example:  Comparing  Interpolating  Functions 


The  following  example  code  creates  a  comma-separated  text  file  that  can  be  used  to  create  the 
graph  shown  in  Fig.  6.  The  plot  relating  to  the  CubeInterpO  function  relies  on  slopes  that  are 
calculated  using  the  CardinalSlope()  function  (except  for  the  endpoints,  where  the  slopes  are  set 
to  zero). 

#include  <cstdio>// . f reopen ( ) ,  printf ( ) ,  stdout 

#include  "y_interp. h"// . yinterp 

int  main(){//<================COMPARISON  BETWEEN  VARIOUS  INTERPOLATION  FUNCTIONS 

freopen("interp.csv"/'w", stdout);// . redirect  output  to  a  file 

double  X[8]={2,7,13,19,22,23,28,37},Y[8]={2,6,6,2,9,5,4,7}; 
double  m[8]={0}; /*<-*/for(int  i=l;i<7;++i) 
m[i]=ylnterp: :CardinalSlope(X+i,Y+i,0. ); 
printf ("#x,y,x,yl,y2,y3\n"); 
for(int  k=0,n=20000;l«n;++k){ 
double  x=50/(n-l. )*k; 

int  i=ylnterp:  :BinarySearch(X+l,X+8-l,x)-X;// . note  0>=i<=6 

double  yl=ylnterp: :NNInterp(X+i,Y+i,x); 

double  y2=ylnterp: : LinInterp(X+i,Y+i,x); 

double  y3=ylnterp: :CubeInterp(X+i, Y+i,x,m[ i] ,m[i+l] ) ; 

k<8?printf("%f,%f/',X[k],Y[k]): printf 

printf ("%f,%f,%f,%f\n",x,yl,y2,y3);} 

^  j  j  E  N  AUT  ^^GI^AI  L  •  LAS  T  '**'U  PDAT  EO'^21DLJL20 


Fig.  6  Comparison  of  the  NNInterpO,  LinInterpO,  and  CubeInterpO  functions 


10.  Example:  Performing  Periodic  Interpolations 


The  following  example  code  begins  by  setting  up  a  set  of  points  that  are  used  to  represent  a 
periodic  function  with  a  period  equal  to  10.0.  Next,  the  NNInterpO,  LinInterpO,  and  CubInterpO 
functions  are  used  to  create  values  for  a  comma-separated  text  file.  The  graph  shown  in  Fig.  7 
represents  the  contents  of  the  output  file. 
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#include  <cstdio>// . f reopen ( ) ,  printf ( ) ,  stdout 

#include  "y_interp. h"// . yinterp 

int  main(){//<=====COMPARISON  BETWEEN  VARIOUS  INTERPOLATION  FUNCTIONS  (PERIODIC) 

freopen( "interp_periodic . CSV" , "w", stdout) ; // . redirect  output  to  a  file 

double  X[4]={10,14,16,18},Y[4]={6,2,9,5}; 

double  XP[7]={X[3]-10,X[0],X[1],X[2],X[3],X[0]+10,X[1]+10}; 
double  YP[7]={Y[3],Y[0],Y[1],Y[2],Y[3],Y[0],Y[1]}; 
double  m[5];/*<-*/for(int  i=0;i<5;++i) 

m[i]=ylnterp: :CardinalSlope(XP+l+i,YP+l+i,0. ); 
printf ("#x,y,x,yl,y2,y3\n"); 
int  j=0; 

for (double  x=-10, k;x<70;x+= . l,++j ){ 

int  i=ylnterp: : PeriodicSearch(X,X+4, k=x,10. ) -X; 

double  yl=ylnterp: :NNInterp(XP+l+i,YP+l+i,k); 

double  y2=ylnterp: : LinInterp(XP+l+i, YP+l+i, k) ; 

double  y3=ylnterp: :CubeInterp(XP+l+i,YP+l+i, k,m[i] ,m[i+l] ); 

j<4?printf("%f,%f/',X[j],Y[j]):printf("-,-/'); 

printf ("%f,%f,%f,%f\n",x,yl,y2,y3);} 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  PDAT  E  1 D  U  L  2  0 


Fig.  7  NNInterpO,  LinInterpO,  and  CubeInterpO  functions  used  to  create  periodic  plots 


11.  Example:  Determining  Interpolation  Performance 


The  following  example  ean  be  used  to  test  the  performance  of  the  NNInterpO,  LinInterpO,  and 
CubeInterpO  functions. 

The  output  was  generated  by  compiling  the  code  using  Microsoft's  Visual  Studio  C++  2010 
Express  compiler,  with  the  output  set  to  “release”  mode.  Surprisingly,  for  this  scenario,  the 
LinInterpO  and  CubeInterpO  functions  outperform  the  NNInterpO  function. 
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#include  <cstdio>// . printf() 

#include  <ctime>// . clock( ) ,CLOCKS_PER_SEC 

#include  "y_interp. h"// . yinterp 

#include  "y_random. h"// . yRandom 

int  main(){//<====PERFORMANCE  COMPARISON  BETWEEN  VARIOUS  INTERPOLATION  FUNCTIONS 
const  int  N=100000000; 
double  X[2]={0,1},Y[2]={0,10}; 

unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 
double  R[N];/*<-*/for(int  i=0;i<N;++i)R[i]=yRandom: :RandU(I,0,l); 
double  S=0,t=clock( ); 

for (int  i=0;i<N;++i)S+=yInterp: :NNInterp(X,Y,R[i] ); 
printf( "Using  NNInterp( ) : \n  y_avg=%.6f,t_avg=%.2f  ns\n\n",S/N, 
(clock()-t)/CLOCKS_PER_SEC/N*lE9),t=clock(),S=0; 
yRandom: : Initialize(1, 1) ; 

for(int  i=0;i<N;++i)S+=yInterp: : LinInterp(X, Y, R[i] ) ; 
printf( "Using  Linlnterp( ) : \n  y_avg=%.6f ,t_avg=%.2f  ns\n\n",S/N, 
(clock()-t)/CLOCKS_PER_SEC/N*lE9),t=clock(),S=0; 
yRandom: : Initialize(1, 1) ; 

for (int  i=0; i<N;++i)S+=yInterp: :CubeInterp(X, Y, R[i] , 10. , 10. ) ; 
printf( "Using  Cubelnterp( ) : \n  y_avg=%.6f ,t_avg=%. 2f  ns\n\n",S/N, 
(clock()-t)/CLOCKS_PER_SEC/N*lE9); 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  PDAT  E  1 D  U  L  2  0 


OUTPUT: 


Using  NNInterpO  ^ 

y_a vg=4 . 999344, t_avg=4 

52 

ns 

Using  LinInterpO: 

y_avg=4 . 999774, t_avg=2 

50 

ns 

Using  CubeInterpO: 

y_avg=4 . 999709 , t_avg=3 

43 

ns 

12.  Example:  Interpolating  in  Two  Dimensions 


The  NNInterpO,  LinInterpO,  and  CubeInterpO  functions  can  be  used  to  perform  interpolations  in 
2  (and  higher)  dimensions. 

Suppose  that  .  =  z{x^ ,  y  . ) .  To  interpolate  in  2  dimensions,  begin  by  performing  2 
interpolations  in  the  y  direction,  1  at  and  1  at  (shown  in  green  in  Fig.  8).  Define  to 
be  the  interpolation  at  x,  and  z^,  to  be  the  interpolation  at  .  Complete  the  interpolation  by 
interpolating  between  z„  and  z^,  in  the  jc  direction  (shown  in  red  in  Fig.  8). 


19 


Fig.  8  Interpolating  in  2  dimensions 

The  following  example  code  uses  the  yBmp  namespace  to  create  3  images  (presented  in  Fig.  9) 
that  show  the  results  of  using  the  NNInterpO,  LinInterpO,  and  CubeInterpO  functions  to 
interpolate  in  2  dimensions. 


#include  "y_interp. h"// . yinterp, <cmath>{fabs( )} 

#include  "y_bmp.h"// . yBmp 

inline  void  Rainbow(//<========================================RAINBOW  COLOR  MAP 

unsigned  char  C[3],//< - OUTPUT  COLOR  (CALCULATED) 

double  x,//< - VALUE  FOR  WHICH  A  COLOR  WILL  BE  CALCULATED 

double  min, double  max){//< - MINIMUM  AND  MAXIMUM  SCALED  VALUES 

if (x<min){C[0]=C[l]=C[2]=0; /*&*/return; }// . set  too  small  values  to  black 

if (x>max){C[0]=C[l]=C[2]=255; /*&*/return; }// . set  too  large  values  to  white 

x=(l-(x-min)/(max-min))*8;// . remap  x  to  a  range  of  8  to  0 

C[0]=int((3<x&&x<5| |x>7  ?-fabs(x/2-3)+1.5:5<=x&&x<=7?l:0)*255);// . blue 

C[l]=int((l<x&&x<3i i5<x&&x<7?-fabs(x/2-2)+1.5:3<=x&&x<=5?l:0)*255);// _ green 

C[2]=int((  x<li i3<x&&x<5?-fabs(x/2-l)+1.5:l<=x&&x<=3?l:0)*255);// . red 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  PDAT  E  D^l  5  D  U  L  2  0  ^4^^'**'^'**"**' 

int  main(){//<======CREATE  AN  IMAGE  FROM  A  SURFACE  USING  INTERPOLATING  FUNCTIONS 

const  int  m=4, n=4, r=250; // . #  of  x  values,  #  of  y  values,  pixel  scaling 

double  X[4]={0,1,2,3},Y[4]={0,1,2,3}; 
double  Z[m][n]={l, 1,1,0  ,  1,0, 0,1  ,  0,1, 0,1  ,  0,0, 0,0}; 
unsigned  char*I=yBmp: :NewImage(m*r,n*r,255); 
for(int  q=0;q<m*r;++q)for(int  p=0;p<n*r;++p){ 
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double  x=q*(m-l)*l./(m*r-l),y=p*(n-l)*l./(n*r-l); 
int  i=ylnterp: :BinarySearch(X+l,X+m-l,x)-X, 
j=ylnterp: : BinarySearch(Y+l, Y+n-l,y) -Y; 
double  ZI[2]={yInterp: :NNInterp(Y+j,Z[i]+j,y), 
yinterp: : NNInterp(Y+j,Z[i+l]+j,y) }; 
double  z=ylnterp: :NNInterp(X+i,ZI,x); 

Rainbow(yBmp: :GetPixel(I,q, p) , z,  - .2,1.2);  } 
yBmp : : WriteBmpFile ( "nearest_neighbor . bmp" , I ) ; 
for(int  q=0;q<m*r;++q)for(int  p=0;p<n*r;++p){ 

double  x=q*(m-l)*l./(m*r-l),y=p*(n-l)*l./(n*r-l); 
int  i=ylnterp: :BinarySearch(X+l,X+m-l,x)-X, 
j=ylnterp: : BinarySearch(Y+l, Y+n-l,y) -Y; 
double  ZI[2]={yInterp: : LinInterp(Y+j,Z[i]+j,y), 
yinterp: : LinInterp(Y+j,Z[i+l]+j,y)}; 
double  z=ylnterp: :LinInterp(X+i,ZI,x); 

Rainbow(yBmp: :GetPixel(I,q, p) , z, - .2,1.2); } 
yBmp : : WriteBmpFile ( "linear . bmp" , I ) ; 
for(int  q=0;q<m*r;++q)for(int  p=0;p<n*r;++p){ 

double  x=q*(m-l)*l./(m*r-l),y=p*(n-l)*l./(n*r-l); 
int  i=ylnterp: :BinarySearch(X+l,X+m-l,x)-X, 
j=ylnterp: : BinarySearch(Y+l, Y+n-l,y) -Y; 
double  ZI[2]={yInterp: :CubeInterp(Y+j,Z[i]+j,y,0. ,0. ), 
yinterp: :CubeInterp(Y+j,Z[i+l]+j,y,0. ,0. )}; 
double  z=ylnterp: :CubeInterp(X+i,ZI,x,0. ,0. ); 

Rainbow(yBmp: :GetPixel(I,q, p) , z, - .2,1.2); } 
yBmp : : WriteBmpFile ( "cubic . bmp" , I ) ; 

^  j  j  E  N  AUT  ^^GI^AI  L  •  LAS  T  PDAT  E  1 D  U  L  2  0 


13.  Polynomial  Curve  Fitting:  The  PolyFit()  Function 


Suppose  that  some  set  of  n  ordered  pairs  y^,)  represents  a  set  of  measurements,  with 
being  a  value  for  an  independent  variable  and  being  a  value  for  a  dependent  variable. 
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Furthermore,  suppose  that  Eq.  22  represents  a  best-fit  equation  for  the  ordered  pairs,  where  the 
coefficients  ( c, )  are  unknown. 


y{x)  =  Cq  -I-  CjX  +  C2X^  H - h  c^x'  H - h  . 

Eq.  23  can  be  used  to  find  the  coefficients  for  Eq.  22:^ 
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The  PolyEitO  function  uses  Gaussian  elimination  with  backward  substitution  to  solve  Eq.  23. 

13.1  PolyFitQ  Code 

template<class  T>void  PolyFit(//<=========FITS  A  POLYNOMIAL  TO  A  SET  OF  POINTS 


const  T*X,//< - INDEPENDENT-VARIABLE  VALUES  (EACH  X[i]  MUST  BE  UNIQUE) 

const  T*Y,//< - DEPENDENT-VARIABLE  VALUES  (SAME  SIZE  AS  X) 

int  n,//< - NUMBER  OF  ELEMENTS  IN  X  OR  Y 

int  m,//< - NUMBER  OF  COEFFICIENTS  IN  THE  BEST-FIT  POLYNOMIAL 

T*C){//< - STORAGE  FOR  COEFFICIENTS  (SIZE=m)  y=C[0]+C[l]*x+C[2]*x''2+. . . 


T**A=new  T*[m];/*<-*/for(int  i=0, j,k;i<m;++i){// . augmented  matrix 

for(A[i]=new  T[m+1], j=0; j<m+l;++j)A[i] [j]=0; 
for(k=0; k<n;++k)for(A[i] [m]+=pow(X[k] , i)*Y[k] , j=0; j<m;++j ) 
A[i][j]+=pow(X[k],i+j);} 

for (int  i=l, j, k; i<m;++i)for(k=0; k<i;++k)for( j=m; j >=0; --j)// . Gaussian 

A[i][j]-=A[i][k]*A[k][j]/A[k][k];//  elimination 

for (int  i=m-l, j;i>=0; --i)for(C[i]=A[i] [m]/A[i] [i], j=i+l; j<m;++j)//. .backward 
C[i]-=C[j]*A[i] [j]/A[i] [i];//  substitution 

for (int  i=0; i<m; ++i) delete [ ]A[i] ; /*&*/delete [ ]A; 

}//  YAGENAUT@GMAIL.COM - LAST~UPDATED~21DUL2014 - 


13.2  PolyFitO  Template  Classes 

T  T  is  typically  a  floating-point  data  type. 

13.3  PolyFitO  Parameters 

X  X  points  to  an  array  that  is  used  to  store  the  n  values  that  are  specified  in 

Eq.  23,  where  X={  Xq, x, ,  •  •  •,  x^ ,  •  •  ■, x„_j } .  Each  Xj,  value  must  be  unique. 

Y  Y  points  to  an  array  that  is  used  to  store  the  n  values  that  are  specified  in 

Eq.  23,  where  Y={  y^,  y^, ■  •  y, , ■  •  y„_i  } . 

n  n  specifies  n  ,  the  number  of  x^  (or  y^ )  values. 
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m  m  specifies  the  number  of  coefficients  in  the  array  that  is  pointed  to  by  C.  Thus, 

m  =  d  +  \,  where  d  is  the  degree  of  the  fitting  polynomial  (see  Eq.  22). 

C  C  points  to  storage  for  an  array  that  is  used  to  store  the  c,  coefficients 

(C= {  Cq  ,  Cj ,  •  •  ■ ,  c, ,  •  •  ■ ,  } ) .  C  must  point  to  an  array  with  storage  for  at  least  m 

elements.  The  elements  pointed  to  by  C  are  calculated  by  the  PolyFit()  function. 

Note  that  if  m  is  too  large,  or  if  the  values  pointed  to  by  X  are  too  close  together,  then  the 
PolyFitO  function  may  return  coefficients  that  do  not  accurately  describe  a  best-fit  curve.  It  is 
always  best  to  plot  the  best-fit  curve  against  measured  data  to  verify  the  quality  of  the  fit.  It  is 
also  a  good  idea  to  calculate  the  coefficient  of  determination  {R^),  which  is  a  measurement  of 
the  quality  of  the  fit. 
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where 


1  k<n 

y  =  -i:y, 

nk=o 


(25) 


values  are  typically  in  the  interval  [0,1],  with  1  indicating  a  perfect  fit. 


13.4  PolyFitO  Simple  Example 

The  following  example  first  defines  a  polynomial,  then  uses  that  polynomial  to  calculate  a  set  of 
{xk,yk)  ordered  pairs.  The  PolyFitO  function  is  used  to  calculate  the  coefficients  for  a 
polynomial  that  best  fits  the  set  of  ordered  pairs.  The  calculated  coefficients  are  shown  to  be 
identical  (at  least  to  6  decimal  places)  to  the  coefficients  of  the  original  polynomial. 


#include  <cstdio>// . printf() 

#include  "y_interp. h"// . yinterp 

int  main(){//<==============A  SIMPLE  EXAMPLE  FOR  THE  yinterp: :PolyFit()  FUNCTION 

double  a=-308,b=177,c=-33,d=2; 
double  X[10],Y[10]; 
for(int  i=0;i<10;++i){ 

X[i]=i, Y[i]=a+b*i+c*i*i+d*i*i*i; 
printf("%f,%f\n",X[i],Y[i]);} 
double  C[4];/*<-*/yInterp: :PolyFit(X,Y,10,4,C); 
printf ("\n"); 

for (int  i=0; i<4; ++i) printf ( "%f\n",C[i] ); 

^11  E  N  AUT  ^^GI^AI  L  •  LAS  T  '**'U  PDAT  ED'^2rDUL20 
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OUTPUT: 


0 . 000000, -308 . 000000 
1 . 000000, -162 . 000000 
2.000000, -70.000000 
3.000000, -20.000000 
4 . 000000, 0 . 000000 
5 . 000000, 2 . 000000 
6 . 000000, -2 . 000000 
7 . 000000, 0 . 000000 
8 . 000000, 20 . 000000 
9 . 000000, 70 . 000000 

-308 . 000000 
177.000000 
-33 . 000000 
2 . 000000 


14.  Example:  Exponential  Fits  Using  the  PolyFit()  Function 


The  PolyFitO  function  can  be  used  to  fit  nonpolynomial  equations  to  measured  data.  Suppose 
Eq.  26  is  suspected  to  be  a  good  fit  for  some  set  of  data  points 

y  =  .  (26) 


The  PolyFitO  function  can  be  used  to  find  the  parameters  A  and  B  . 

Begin  by  defining  y'  =  ln(y) .  Next,  use  the  PolyFitO  function  to  find  and  Cj  (the  parameters 
for  a  straight  line)  for  the  data  set  k^yO-Thus, 


Solving  for  y , 

Thus, 


y  =  ln(y)  =  Co-l-CiX. 


y  =  e 


A  =  e'"”  and  B  =  c^. 


(27) 


(28) 

(29) 


The  following  example  code  uses  the  PolyFitO  function  to  fit  Eq.  26  to  beginning-of-year  values 
for  the  Dow  Jones  Industrial  Average  for  the  decade  1986-1995.  Values  were  obtained  from  the 
file  “Mji_d.csv,”  which  was  downloaded  from  http://stooq.com.^  The  yIo2  namespace  was  used 
to  read  and  parse  the  data  file.  The  results  are  presented  in  Fig.  10. 
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#include  <cstdio>// . f reopen ( ) ,  printf ( ) ,  stdout 

#include  <vector>// . vector 

#include  "y_interp. h"// . yinterp, <cmath>{exp( ) , pow( )} 

#include  "y_io_2.h"// . ylo 

int  main(){//<======FITTING  AN  EXPONENTIAL  EQUATION  USING  THE  PolyFit()  FUNCTION 

// - read  and  parse  file - 

char*s=yIo2: : ReadTextFile( "^d ji_d . csv" ); 

std: :vector<std: :vector<char*>  >S=yIo2: :Parse2D(s, "/ 

// - INTERPRET  DATES  AND  VALUES - 

std; : vector<double>X,Y; 

for (int  i=2,y; (y=atoi(S[i] [0] ) )<1996;++i) 

if (y ! =atoi(S[i-l] [0] )&&y>1984)X.push_back(y) , Y. push_back(atof (S[i] [6] ) ) ; 
delete[ ] s; 

// - PERFORM  EXPONENTIAL  FIT - 

std : : vectorc double >Yp(X. size ( ) ) ; 

for (int  i=0,m=X. size( ) ; i<m;++i)Yp[i]=log(Y[i] ) ; 

double*a=&X[0] , *b=&Yp[0] ; 

double  C[2];/*<-*/yInterp: :PolyFit(a,b,X.size(),2,C); 

double  y_bar=0; /*<-*/for(int  k=0; k<10;++k)y_bar+=Y[k] ; /*&*/y_bar/=10; 

double  SS_res=0,SS_tot=0; /*<-*/for(int  k=0; k<10;++k) 

SS_res+=pow(Y[k] -exp(C[0]+C[l]*X[k] ) ,2. ),SS_tot+=pow(Y[k] -y_bar,2. ); 
printf ( "A=%E\nB=%f\nR^2=%f\n",exp(C[0] ) ,C[1] , l-SS_res/SS_tot) ; 
f reopen ( "exponential_f it .csv" , "w", stdout) ; 
for(int  i=l,n=X.size();i<n;++i) 

printf ("%f,%f,%f\n",X[i],Y[i],exp(C[0]+C[l]*X[i])); 

^11  E  N  AUT  ^)GI^AI  L  •  LAS  T  PDAT  E  D^2 1 D  U  L  2  0 


Fig.  10  Best-fit  line  calculated  by  the  PolyFit()  function 
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15.  Code  Summary 


The  following  summary  sheet  presents  the  yinterp  namespace,  which  contains  the 
BinarySearchO,  PeriodicSearch(),  NNInterpO,  LinInterpO,  CubeInterpO,  CardinalSlope(),  and 
PolyFitO  functions. 
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yinterp  Summary 


y_interp.h 


#ifndef  Y_INTERP_GUARD//  See  Yager,  R.D.  "Basic  Searching,  Interpolating,  and 
#define  Y_INTERP_GUARD//  Curve  Fitting  Algorithms  in  C++"  (ARL-TN-XXX) 

#include  <cmath>// . fmod()jpow() 

namespace 

templatecclass  T>T*BinarySiffW(//<======FIND  THE  POINTER  c  |  *c  <=  k  <  *(c+l) 

T*a,T*b,//<--ARRAY  START  &  END  POINTS  (ARRAY  MUST  BE  SORTED  IN  INC.  ORDER) 

T  k){//< . KEY 

if (k<*a)return  a-1;// . note  that  a-1  may  point  to  an  invalid  address 

for(T*c; k<*--b; k>*c?a=c :b=c+l)c=a+(b-a)/2j /*->*/ return  b; 

YAGENAUT@GMAIL.COM - LAST~UPDATED~213UL2014 - 

templatecclass  T>T»priodicSearcH(//<===*===FOR  ARRAYS  WITH  PERIODIC  VARIABLES 
T*a,T*b,//<--STARTING  &  ENDING  POINTS  (ARRAY  MUST  BE  SORTED  IN  INC.  ORDER) 

T&k,//< . KEY  (WILL  BE  SET  TO  k'=k+np) 

T  p){//< . PERIOD  (p>0) 

return  Binary Search ( a, b,k=f mod (k-*a,p)+*a+(k-*a<0?p:0))] 

YAGENAUT@GMAIL.COM - LAST-UPDATED~2i:UL2014 - 

template<class  T>T  NNInterp(//<============«====NEAREST-NEIGHBOR  INTERPOLATOR 

const  T*X,//< - BRACKETING  X  VALUES  (*X  AND  X[l]  MUST  BE  VALID) 

const  T*Y,//< . BRACKETING  Y  VALUES  (*Y  AND  Y[l]  MUST  BE  VALID) 

T  x){//< . VALUE  TO  INTERPOLATE  AT  (TYPICALLY,  *X  <=  X  <=  X[l]) 

return  x>(»X+X[l])/2?Y[l] :*Y; 

YAGENAUT@GMAIL.COM— - - - — - - LAST~UPDATED'-213UL2014 - — 

templatecclass  T>T  Linlnterp(//<===========«=*====**===«=LINEAR  INTERPOLATOR 

const  T*X,//< - BRACKETING  X  VALUES  (*X  AND  X[l]  MUST  BE  VALID) 

const  T*Y,//< . BRACKETING  Y  VALUES  (*Y  AND  Y[l]  MUST  BE  VALID) 

T  x){//< . VALUE  TO  INTERPOLATE  AT  (TYPICALLY,  *X  <=  X  <=  X[l]) 

return*Y+(Y[l]-*Y)*(x-*X)/(X[l]-*X); 

}// - YAGENAUT@GMAIL.COM - - - — - — - - LAST~UPDATED'-213UL2014 - — 

templatecclass  T>T  tubeInter^(//<=======*==CUBIC  (HERMITE  SPLINE)  INTERPOLATOR 

const  T*X, //<--- - BRACKETING  X  VALUES  (*X  AND  X[l]  MUST  BE  VALID) 

const  T*Y,//< . BRACKETING  Y  VALUES  (*Y  AND  Y[l]  MUST  BE  VALID) 

T  x,//< . VALUE  TO  INTERPOLATE  AT  (TYPICALLY,  *X  <=  X  <=  X[l]) 

T  m0,T  ml){//< . SLOPES  AT  *X  AND  X[l] 

T  t=(x-*X)/(X[l]-*X); 

return*Y+(x-*X)*(((m0+ml)*t-(2*m0+ml))*t+m0)+(Y[l]-*Y)»(3-2*t)*t*tj 

}// - YAGENAUT@GMAIL.COM~ - - - - - — — — LAST~UPDATED~2i:UL2014~~~~~~ 

template<class  T>T  CardinalSlope(//<====®===CALCULATES  SLOPES  FOR  CubeInterpO 

const  T*X,//< - BRACKETING  X  VALUES  (X[-l]  AND  X[l]  MUST  BE  VALID) 

const  T*Y,//< . BRACKETING  Y  VALUES  (Y[-l]  AND  Y[l]  MUST  BE  VALID) 

T  t){//< . TENSION  PARAMETER 

return(l-t)*(Y[l]-Y[-l])/(X[l]-X[-l])j 

}// - YAGENAUT@GMAIL.COM~ - - - - - ~~~~~~~~~LAST~UPDATED~213UL2014~~ - 

templatecclass  T>void  PolyFit(//<=========FITS  A  POLYNOMIAL  TO  A  SET  OF  POINTS 

const  T*X,//< . INDEPENDENT -VARIABLE  VALUES  (EACH  X[i]  MUST  BE  UNIQUE) 

const  T*Y,//< . DEPENDENT-VARIABLE  VALUES  (SAME  SIZE  AS  X) 

int  n,//< . NUMBER  OF  ELEMENTS  IN  X  OR  Y 

int  m,//< . NUMBER  OF  COEFFICIENTS  IN  THE  BEST-FIT  POLYNOMIAL 

T*C){//< . STORAGE  FOR  COEFFICIENTS  (SIZE=m)  y=C[0]+C[l]*x+C[2]*x''2+. . . 

T**A=new  T*[m];/*<-*/for(int  i=0, jjk;i<mj++i){// . augmented  matrix 

for(A[i]=new  T[m+1] ,j=0j j<m+lj++j)A[i] [j ]=0j 
for(k=0;k<nj++k)for(A[i] [m]+=pow(X[k]ji)*Y[k], j=0; j<m;++j) 
A[i][j]+=pow(X[k],i+j);} 

for (int  i=l, jjk;i<mj++i)for(k=0;k<i;++k)for(j=mj j>=0;--j)// . Gaussian 

A[i][j]-=A[i][k]*A[k][j]/A[k][k];//  elimination 

for (int  i=m-lj j;i>=0; --i)for(C[i]=A[i] [m]/A[i] [i] , ]=i+lj ] <m; ++])// . . backward 
C[i]-=C[j]*A[i] [ j ]/A[i] [i] ; //  substitution 

for (int  i=0;i<m;++i)delete[]A[i] j/*&*/delete[]A; 

}//.w....,~YAGENAUT@GMAIL.COM~~ - - - - - - - ~~~LAST~UPDATED~213UL2014~ - ~~ 


#endif 


BinarySearchO  Performance 


#include  <algorithm>// . sort() 

#include  <cstdio>// . printf() 

#include  <ctime>// . clock( ),CLOCKS_PER_SEC 

#include  "y_interp.h”// . yinterp 

#include  "y_randoffl.h”// . yRandom 

templatecclass  T>T*LinearSearch(//<========FIND  THE  POINTER  c  |  *c  <=  k  <  *(c+l) 

T*a,T*b,//< - ARRAY  START  &  END  POINTS  (ARRAY  MUST  BE  SORTED  IN  INC.  ORDER) 

T  k){//< . KEY 

if(k<*a)return  a-1;// . note  that  a-1  may  point  to  an  invalid  address 

while(k<*--b);/*->*/return  b; 

}//.w - YAGENAUT@GMAIL.COM~~ - - - - - ~~~~~~~~LAST~UPDATED~2i:UL2014~ - 

int  main(){//<==*«====PERFORMANCE  TEST  FOR  THE  yinterp: : BinarySearch( )  FUNCTION 
const  int  N=16384,M=10000000; // . . .max  #  of  elements  in  array,  #  of  repetitions 
unsigned  I[625];/*<-*/yRandom: :Initialize(I,l);//. .. .state  of  Mersenne  twister 
double*X=new  double[N];/*<-*/for(int  i=0;i<N;++i)X[i]=yRandom: :RandU(Ij0,N); 
std: :sort(XjX+N);// . X  is  a  sorted,  unchanging  list  of  numbers  to  search 


array 


search  | 
time  (s)  I 

-I- 


index 

average 


binary 


search  | 
time  (s)  I 
- 1 


search\n"); 
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index\n"); 
average\n"); 
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}//• 


printf("  I  linear  search 

printf ( " 
printf ( " 
printf ( " 
printf ( " 

for(int  m=2;m<16385;m*=2){ 

printf ( "%8d  | ",m), yRandom: :Initialize(Ij 1); 

double  y=0,z=0,t=clock(); 

for (int  i=0 ;i<M;++i)y+= Linear Search (X,X+m, yRandom: :RandU(Ij0,m) ) -X; 
printf ("%8.3f  |%12.6f  | ", (clock( ) -t)/CLOCKS_PER_SEC,y/M) ,t=clock( ); 

yRandom : : Initialize(Ij 1); 

for (int  i=0;i<M;++i)z+=yInterp: : BinarySearch(X,X+m, yRandom: : RandU(I,0jm) ) -X; 
printf("%8.3f  |%12.6f\n",(clock()-t)/CL0CKS_PER_SEC,z/M);} 

- — YAGENAUT@GMAIL.COM~~ - - - ~~~~LAST~UPDATED~2i:UL2014~ - 


Linear-Search  Performance  -  0{n) 


Binary-Search  Performance  -  0(log(n)) 


— 

BinarySearch()/PeriodicSearch( )  Example 


#include  <cstdio>// . freopen( ), printf (), stdout 

#include  "y_interp.h"// . yinterp, <cmath>{sin()} 

int  main( ){//<====**=====®==COMPARISON  BETWEEN  BinarySearchO  &  PeriodicSearch() 
f reopen ("sine. CSV", "w", stdout); 

double  X[20];/*<-*/for(int  i=0;i<20;++i)X[i]=5*i/20.+5; 
double  Y[20];/*<-*/for(int  i=0;i<20;++i)Y[i]=sin(2*3.14159*X[i]/5)+l; 
printf("x,y  -  BinarySearchO), y  ■  PeriodicSearchO\n''); 
for(double  x=-5,k;x<15;x+=.l){ 

printf ( "%f ,%f , ”,x,Y [yinterp: : BinarySearch(X+ljX+20,x) -X] ) ; 
printf ( "%f\n",Y [yinterp: : PeriodicSearch(X,X+20, k=x,5. ) -X] ); } 

}//— - YAGENAUT@GMAIL.COM~~~~~~~~~ - - - ~~~~~~~LAST~UPDATED~2i:UL2014~~ - - 
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Interpolation  Comparison 


#include  <cstdio>// . f reopen (), printf (), stdout 

#include  "y_interp. h"// . yinterp 

int  main(){//<*====«==«=====COMPARISON  BETWEEN  VARIOUS  INTERPOLATION  FUNCTIONS 

freopen(”interp. CSV"  J  ”w",  stdout); // . redirect  output  to  a  file 

double  X[8]={2,7,13,19,22,23,28,37},Y[8]={2,6,6,2,9,5,4,7}; 
double  m[8]={0};/*<-*/for(int  i=l;i<7;++i) 
m[i]=ylnterp: :CardinalSlope(X+i,Y+ij0. ) ; 
printf("#XjyjX,yljy2,y3\n"); 
for(int  k=0,n=20000;k<n;++k){ 
double  x=50/(n-l. )*k; 

int  i=ylnterp:  :BinarySearch(X+ljX+8-l,x)-X;// . note  0>=i<=6 

double  yl=ylnterp: :NNInterp(X+i,Y+ijX); 

double  y2=ylnterp: : LlnInterp(X+i,Y+i,x); 

double  y3=ylnterp: :CubeInterp(X+ijY+i,x,m[i],m[i+l]); 

k<8?printf("%f,%f,",X[k],Y[k]): printf ("-,-,"); 

printf("%f,%f.%f,%f\n",x,yl,y2,y3);} 


PolyFitO  Example 


#include  <cstdio>// . freopenO, printf (), stdout 

#include  <vector>// . vector 

#include  "y_interp. h"// . yinterp, <cmath>{exp() , pow()} 

#include  "y_io_2.h”// . ylo 

int  main(){//<======FITTING  AN  EXPONENTIAL  EQUATION  USING  THE  PolyFit()  FUNCTION 

// . READ  AND  PARSE  FILE . 

char*s=yIo2 : :  ReadTextFile("''dji_d .  csv" ); 

std: :vector<std: :vector<char*>  >S=yIo2: : Parse2D(s, "/  \t-,"); 

// . INTERPRET  DATES  AND  VALUES . 

std: :vector<double>XjY; 

for (int  i=2,y; (y=atoi(S[i] [0] ) )<1996;++i) 

if (y ! =atoi(S[i-l] [0] )&&y>1984)X.push_back(y) ,Y.push_back(atof (S[i] [6] ) ); 
delete[]s; 

// . PERFORM  EXPONENTIAL  FIT . 

std : : vector<double>Yp(X.size( )) ; 

for (int  i=0jm=X.size();i<m;++i)Yp[i]=log(Y[i]); 

double*a=&X[0] , *b=&Yp[0] ; 

double  C[2];/*<-*/yInterp: :PolyFit(a,bjX.size()j2,C); 

double  y_bar=0;/*<-*/for(int  k=0;k<10;++k)y_bar+=Y[k];/*&*/y_bar/=10; 

double  SS_res=0,SS_tot=0;/*<-*/for(int  k=0;k<10;++k) 

SS_res+=pow(Y[k] -exp(C[0]+C[l]*X[k] ) J  2. ) ,SS_tot+=pow(Y[k] -y_bar,2. ); 
printf ("A=%E\nB=%f\nR^2=%f\n",exp(C[0]),C[l],l-SS_res/SS_tot); 
freopen("exponential_fit.csv",”w"j stdout); 
for(int  i=ljn=X.size();i<n;++i) 

printf ("%f,%f,%f\n",X[i],Y[i],exp(C[0]+C[l]*X[i])); 

}// - YAGENAUT@GMAIL.COM~~ - - - LAST~UPDATED~2i:UL2014~~~ - 
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